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In previous work we have shown that the Density Matrix Renormalization Group (DMRG) en- 
ables near-exact calculations in active spaces much larger than are possible with traditional Complete 
Active Space algorithms. Here, we implement orbital optimisation with the Density Matrix Renor- 
malization Group to further allow the self-consistent improvement of the active orbitals, as is done in 
the Complete Active Space Self-Consistent Field (CASSCF) method. We use our resulting DMRG- 
CASSCF method to study the low-lying excited states of the all-trans polyenes up to C24H26 as well 
as /3-carotene, correlating with near-exact accuracy the optimised complete 7r-valence space with up 
to 24 active electrons and orbitals, and analyse our results in the light of the recent discovery from 
Resonance Raman experiments of new optically dark states in the spectrum. 
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I. INTRODUCTION 

The Density Matrix Renormalization Group (DMRG) 
is an eleetronic structure method that has recently been 
applied to ab-initio quantum chemistry. The method 
originated in the condensed matter community with the 
pioneering work of White [1, 2]. Although the earli- 
est quantum chemistry implementations are only a few 
years old, the DMRG has already enabled the solution 
of many problems that would be intractable with any 
other method [3-7]. For example, we have shown that 
the DMRG can obtain near-exact solutions to multiref- 
erence problems with active spaces much larger than are 
possible with traditional active space techniques. Such 
problems have ranged from molecular potential energy 
curves [8, 9], to the ground and excited states of large 
conjugated polymers [10-12], to metal-insulator transi- 
tions in hydrogen chains [10]. In each of these cases, 
we obtained DMRG energies within 0.001-0. lm£^/i of the 
(estimated) exact Full Configuration Interaction (FCI) 
energies in the active space, but for active spaces that, in 
some problems, have been as large as 100 active electrons 
in 100 orbitals [10]. The development of the DMRG in 
quantum chemistry has proceeded through the efforts of 
several groups, and we mention here the work of White 
et al. [3, 13, 14], Mitrushenkov et al. [4, 15, 16], our 
contributions [5, 8-12, 17, 18], the work of Legeza, Hess 
et al [6, 19-21], the work of Reiher et al. [7, 22-24], and 
most recently the work of Zgid and Nooijen [25]. Also 
related, but too numerous to cite in full here, are earlier 
developments of the method for semi-empirical Hamil- 
tonians; some representative contributions are those in 
Refs. [26-32]. 

At the heart of the DMRG is a wavefunction ansatz 
and the DMRG "algorithm" is simply an efficient vari- 
ational optimisation procedure for this ansatz. Unlike 
most wavefunctions in quantum chemistry, the DMRG 



wavefunction is not parametrised by excitations from an 
underlying reference state. Rather, it is built directly 
from local variational objects (which we shall later call 
site functions) which are associated with the active or- 
bitals in the system, and which describe how the orbitals 
are correlated with each other. Each site function is 
characterised by a rank M that measures the number 
of variational parameters, and as this rank increases the 
ansatz becomes exact. For an incomplete rank M, cor- 
relations between orbitals that are widely separated in 
the ansatz are truncated. Thus the DMRG is a natu- 
rally local theory, but, since the ansatz is not constructed 
from a reference, it is a local multireference theory. This 
may be seen as the basic reason why the DMRG can de- 
scribe very large multireference problems so easily. We 
should note that the structure of the DMRG wavefunc- 
tion means that it is a local theory only in the number of 
correlating orbitals along one of the physical dimensions 
of the problem. However, generalisations of the ansatz 
to a local theory along all physical dimensions are now 
known, and are under active development [33-38]. 

In most applications of the DMRG to quantum chem- 
istry so far, the active space of interest has been easy 
to identify, i.e. there is a good core-valence and valence- 
Rydberg separation, either for energetic or symmetry rea- 
sons, allowing the DMRG to be used with such an active 
space as a direct substitute for Complete Active Space 
Configuration Interaction (CASCI). In general, however, 
we cannot always identify the active orbitals in a sim- 
ple way, and thus there is a need for an orbital opti- 
mised DMRG, where the active space is determined self- 
consistently by energy minimisation, in much the same 
way as in the Complete Active Space Self-Consistent 
Field (CASSCF) method [39, 40]. The purpose of the 
current work is to describe how this may be done. The 
resulting orbital optimised DMRG we shall refer to as 
the DMRG-CASSCF method. 
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While the general idea of orbital optimisation is 
straightforward, in practice an efficient implementation 
must be tailored to the underlying many-body wavefunc- 
tion ansatz. In Sec. II we describe such an algorithm for 
the DMRG wavefunction. We start with an overview of 
orbital optimisation in Sec. II A that recalls how the pro- 
cedure may naturally be divided into two tasks, the eval- 
uation of the one- and two-particle density matrices, and 
the orbital rotation and integral transformation steps. In 
Sec. II B we present an efficient method to evaluate the 
one- and two-particle density matrices in the DMRG. Our 
current implementation benefits from the observation of 
Zgid and Nooijen that the one-site DMRG algorithm is 
more suitable than the two-site DMRG algorithm for this 
purpose. To facilitate the large-scale calculations for our 
applications to long polyenes and /3-carotene in this work, 
we have fully parallelised not only the evaluation of the 
reduced density matrices in the DMRG, but also the or- 
bital rotation and integral transformation steps. These 
implementation aspects are discussed in Sec. II C. Fi- 
nally, the complete DMRG-CASSCF macroiteration is 
summarised in Sec. II D. 

In Sec. Ill we apply the DMRG-CASSCF method to 
the problem of the low-lying excitations in polyenes and 
/3-carotene. The conjugated 7r-system in the polyenes 
and substituted species such as /3-carotene gives rise to 
an unusual excitation spectrum, with "dark" electronic 
states lying beneath the optically allowed HOMO-LUMO 
transition. The electronic structure of these low-lying 
states lies at the heart of energy transport in systems 
ranging from conjugated organic semiconductors to the 
biological centres of light-harvesting and vision. While 
the relevant active space on these systems clearly consists 
of the conjugated 7r-valence orbitals. to the best of our 
knowledge previous calculations on these systems have 
not correlated complete tt- valence spaces with more than 
5 double bonds (corresponding to a (10,10) complete ac- 
tive space [41, 42]). In the current study we use our 
DMRG-CASSCF method to perform calculations corre- 
lating the complete 7r-valence space in polyenes up to 
C24H26 (with 12 conjugated bonds) and /3-carotene (with 
11 conjugated bonds), and analyse our results in relation 
to recent Resonance Raman measurements, which have 
detected previously unidentified "dark" states in the low- 
lying spectrum. 



II. THEORY 

A. Overview of orbital optimisation 

We begin with some general remarks on orbital optimi- 
sation in ab-initio quantum chemistry. Starting from the 

elcKitronic Hamiltonian, specified by the one- and two- 
electron integral matrix elements Uj and Vijki 

if = ^ Ujolaj + ^ v,jkiala]jakai (1) 

ij ijkl 



an ab-initio quantum chemical method provides a wave- 
function that approximates a target eigenstate of H. 
From \E' we define the one- and two-particle density ma- 
trix elements jij,^ijki 

7ii = (*l«I«,l*) (2) 
lijki = {^\ala]akai\^) (3) 

and the energy expectation value (^|ii|^') can be written 
as 

-E' = X] + X] ^ijknijkl (4) 

ij ijkl 

Orbital rotation corresponds to a unitary transforma- 
tion of the wavefunction effected by an operator e^, 
where A has the single-particle operator form 

A = J2^ij4aj (5) 

ij 

and Aij = —A*,-. After orbital rotation, the transformed 
wavefunction ^ and energy E are 

E={^e-^\H\e^^) (6) 

But one can also consider the unitary operator to act on 
the Hamiltonian rather than the wavefunction, and from 
this equivalent point of view, we have a transformed H 
and energy expression 

H = e-^He^ 

E=mH\<i>) (7) 

The transformed Hamiltonian H has the same form as 
the original Hamiltonian (1) but with modified integrals 
tij and Vijki that reflect the rotated orbitals 

iij = U^^l Ujji ti'ji 

i'j' 

Vijki = X '^ii'^jV^'sfe'^«''"j'j'fe'(' (^) 

i'j'k'l' 

where U is the coefficient matrix . Thus we can rewrite 
the energy after orbital rotation in terms of the original 
one- and two-particle density matrices and the modified 
integrals 

= X] *y Tii + X] "^ijkHijkl (9) 
ij ijkl 

We include this elementary discussion because it leads 
directly to the following familar procedure to optimise 
the orbitals in an ab-initio wavefunction: 

1. Prom the ab-initio method obtain ^ correspond- 
ing to the given H and form the density matrices 

^ij 5 ^ijkl ' 
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2. Determine an orbital rotation step eA, and form 
the new Hamiltonian H = e~^He^ from the trans- 
formed integrals. 

3. Goto 1. and loop until convergence in ^. 

Note that in the above, the orbital degrees of freedom and 
the other ansatz degrees of freedom in \1/ are alternately 
optimised in steps (1), (2). While more sophisticated ap- 
proaches which couple orbital rotations with changes in 
the other ansatz degrees of freedom can be envisaged (as 
are employed in multi-coiifigurational self-consistent field 
methods [43, 44]), we shall adopt the above simple strat- 
egy to optimise the orbitals in the DMRG wavefunction. 
The conceptual task is then twofold. Firstly, how do we 
calculate the one- and two-particle density matrices in 
the DMRG? And secondly, what method should we use 
to select our orbital rotation steps and to construct the 
transformed Hamiltonian? 



B. Evaluation of the one- and two-particle density 
matrices in the DMRG 

While the algorithm to calculate the one- and two- 
particle density matrices could, in principle, be described 
entirely in the traditional Renormalization Group lan- 
guage of the DMRG, we believe that it is beneficial to 
understand the method in a more modern language which 
focuses on the structure of the DMRG wavefunction. 
Thus we begin with a brief review of the general prop- 
erties of the DMRG wavefunction before proceeding to 
the method of reduced density matrix evaluation. For an 
expanded introduction to the wavefunction perspective 
in DMRG, we refer the reader to our introductory article 
Ref. [45] as well as other recent reviews in the field [46]. 

1. The DMRG wavefunction 

The DMRG algorithm corresponds to a variational 
minimisation of the energy within the space of a wave- 
function ansatz. To specify this ansatz we first define 
an ordering of the orbitals thereby mapping them onto 
sites on a one-dimensional lattice. Then, the "one-site" 
DMRG ansatz is given by 

I^dmrg) = Yl CCi. Cis • • • C-1 1"i"2n3 ... rife 

nin2n3...nk 
iiigia-.-ife-i 

(10) 

where \ni . . .Uk) denotes a Slater determinant in occu- 
pation number form, i.e. rij is the occupation of orbital 
i, and the total number of orbitals is k. The ip "site 
functions" are 3-indcx quantities and are the variational 
parameters of the wavefunction. The dimension of each 
ni . . . rifc index is 4, corresponding to the 4 occupancies 
for each orbital |— ), [</)"), \4'^), \4>"(f>^), while the dimen- 
sion of each auxiliary index ii . . .ik-i is some specified 



size M, thus making each site function a tensor of di- 
mension 4 X M X M, except for the first and last, which 
only have two indices and are of dimension 4 x M. As 
M increases, the wavefunction ansatz becomes increas- 
ingly exact. If we interpret a site function with indices 
np,ip-i,ip as a matrix array ip"'' where ip-i,ip are the 
matrix indices and rip is the third, array, index, then the 
ansatz is written compactly as a matrix product state 

|*dmrg)= Y1 V'"^V'"=V'"'---V'"Mnin2n3...nfc) 

nin2n3...nk 

(11) 

Because of this matrix product structure, the DMRG 
ansatz is also known as the matrix product state (MPS) 
[47-49]. 

Now the above form of the DMRG ansatz is invari- 
ant to transformations of the site functions of the form 
(^"p ^ ^"pU, ■0"''+' W'ip-^p+^) and thus it is useful 
to define a canonical form of the DMRG wavefunction 
that eliminates this freedom. In practice, this canoni- 
cal representation is used in all DMRG calculations, and 
it is also the representation in which the link between 
the DMRG wavefunction and the traditional Renormal- 
ization Group language is most direct. In essence, the 
canonical form of the wavefunction at a given site corre- 
sponds to the familiar expression for the DMRG wave- 
function where it is expanded in the product basis of the 
left and right blocks separated by the site [3, 45]. 

To obtain the canonical form, we choose a specific site, 
say p, around which to canonicalise. Then the site p 
canonical form is given as 

1^')= J2 i"' •.•-Z^"''-'C'"''-R"''+' •.•-R"'|ni... rip... rife) 

ni...np...nk 

(12) 

We label the site functions to the left of p by i, and those 
to the right by R. The degeneracy (invariance to trans- 
formation) of the original ansatz (10) mentioned above 
is lifted by requiring the L and R site functions to be 
orthogonal projection matrices in the following sense 

Y,L:,JLlJ,=5in" (13) 

Y,K'rK''r = ^r'r" (14) 

rriq 

i.e. by grouping together the Iriq indices to form the row 
index of a 4M x M matrix, each L site function is orthog- 
onal with respect to its M columns, while by grouping 
together the rrig indices to form the column index of a 
M X 4M matrix, each R site function is orthogonal with 
respect to its M rows. 

The link between the canonical form and the orig- 
inal RG formulation appears when we combine the 
L site functions L"^^ . . . L^^-^ with the basis states 
[rii . . . rip-i), and the R site functions i2"p+i . . . i?"*" with 
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the basis states |np+i . . . nt), to define renormalised left 
and right many body spaces {rp+i} 



i^p-i>= E ^r---ir;_".UK---«p-i) (15) 

ni...np_i 

ll...lp-2 

\r,+ ,)= E Rr::^r,+.---Rr^K+l.-.nk) (16) 
np+i...nfe 
rp+2...rfe 

Since the dimension of the left basis in Eq. (15) is M 
(i.e. the dimension of the auxiliary index Ip-i) and sim- 
ilarly for the right basis, the site functions L^^ . . . L^p-'^ 
and . . . ii"** define a projective transformation or 

renormalization from the many-body spaces {ni} (g) . . . (g) 
{up-i} and {rip+i} (gi . . . (g) {nfc} to the left and right 
spaces, {Ip-i} , {rp+i} , respectively. Then, in the renor- 
malised representation, C"'' gives the coefi&cients of 
expansion of the wavcfunction I^P), i.e. 



I*) = 



E 

ip — ]^ 



(17) 



This is just the RG expression for the one-site DMRG 
wavefunction, in the product space of a renormalised 
left "block" , a site p, and a renormalised right "block" . 
Thus in the usual DMRG language, the site p canonical 
form corresponds to the DMRG wavefunction in the basis 
associated with the block configuration 



•p-i-i 



A one-site DMRG wavefunction expressed in the 
canonical form of a given site p can always be expressed 
in the canonical form for any other site (or using the tra- 
ditional DMRG language, the DMRG wavefunction for 
a given one-site block configuration can always be ex- 
pressed in the basis of any other one-site block configu- 
ration along a sweep). Since we are simply re-expressing 
the same wavefunction in a difi^erent basis, the coefficients 
C and site- functions L,Rat difl[erent sites are related. To 
see the link explicitly, we compare the canonical forms at 
adjacent sites p, p+1 



|*)= ...L"''-iC"''H"''+iR"''+=^ ...i2"'=|ni...np...nfe) (18) 

ni...np...nk 

= ^ ...L"f-iL"fC"''+ii?"»'+^ ...i«"'=|ni... rip... rife). (19) 



ni...np...nk 



which yields the relation 



(20) 



Now say we are given C"pi?"p+^ from the site p canoni- 
cal form, and we wish to determine i"pC"p+'^ for the site 
p+1 canonical form, where satisfies the orthogonal- 
ity conditions (13). We can obtain such a L"" solution of 
(20) together with C"p+^ from the singular value decom- 
position (SVD) of C"f , viewed as the 4M x M matrix 
with row indices Ip-iUp, column indices rp+i and M sin- 
gular values cTjp , 

C'C-i.'-p+i = E-^r/_i,ipCrip"^ip'-p+i' (21) 

= E <^i.vi.rp,.Rr::u.. (22) 



The above transformation between canonical forms at 
adjacent sites corresponds directly to the transformation 
between block configurations during the sweep algorithm 
in the DMRG. In particular, Eq. (21) corresponds to 
the determination of the basis of the renormalised block 
•i . . . "p+i I from the density matrix eigenvectors of the 
superblock •! . . . I "p+i, while Eq. (22) corresponds 
to the wavefunction transformation used to generate the 
guess at a given block configuration from that at the 
previous configuration. We note in passing that an exact 
transformation between canonical forms at diS^erent sites 
is only possible with the one-site DMRG ansatz. Most 
DMRG calculations use the two-site DMRG ansatz with 



the block configuration 
and a corresponding canonical form at site p 



'p-i 



•p+2 ■ 



'■p+1 
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|*>= ^ i"i...L"^-iC"''"''+ii?"''+^...i«"*'|ni...npnp+i...nfe) 

ni...np...nk 

Yl C'C-"^rp+2l^P-l'^P'^P+l''P+2) 
/p_ 1 npTip-i- 1 rp-|-2 



(23) 
(24) 



Unlike in the one-site ansatz, the coefficient matrix 
Qnpnp+i jjg^g different shape from the L and R site func- 
tions and has AM (as opposed to M in the one-site case) 
singular values. Thus it can only be approximately rep- 
resented by the sum over M singular values in Eq. (22), 
and the resulting truncation corresponds to "discarding 
states", in the DMRG algorithm. The primary bene- 
fit of the two-site DMRG ansatz is greater robustness of 
convergence in the DMRG sweeps but for the purposes of 
orbital optimisation, the one-site DMRG ansatz provides 
a single consistent DMRG wavefunction in all canonical 
forms and block configurations and is to be preferred. 



2. Reduced density matrix evaluation 

Our task now is, given a DMRG wavefunction writ- 
ten explicitly as (12) or cquivalcntly in the renormalised 
expansion (17), to find an efficient algorithm to evalu- 
ate the one- and two-particle density matrices. Prom 
the renormalised form we see that we will need matrix 
representations of operators in each of the three spaces 
{Zp_i}, {rip}, {rp-i-i}, i.e. matrix elements {P~^\0\P~^ ), 
{nP\0\nP') , {rP^^\0\rP^^ ). Matrix representations in the 
left and right spaces are in general of dimension M x M, 
since there are M left and right states. While the di- 
rect evaluation of the one-particle density matrix would 
require operator representations and thus 0{M^k^) 
storage (presenting no particular difficulties as the mem- 
ory requirement for the usual DMRG algorithm is also 
O(Af^fc^)) the two-particle density matrix would require 
O(M^fc^) storage which is prohibitively expensive. (It 
might appear that when solving the Schrodinger equa- 
tion, the action H\^) would also involve 0{k^) operators 
and 0{APk^) storage. However, there we do not need the 
action of the operators ala^jaka,i individually, but only 

the total ^ijkiVijkio\a}jakai^ so wc can form intermedi- 
ates where operators are precontracted with two-electron 
integrals to save memory, and the efficient arrangement 
of such intermediates lies at the heart of the quantum 
chemical DMRG algorithm). 

The way forward is to observe that we are not tied to 
using a single canonical form/block configuration for the 
DMRG wavefunction, but rather, can evaluate a den- 
sity matrix element ^ijki at any canonical form/block 
configuration that is convenient. As we have described 
above, a given DMRG wavefunction can be expressed in 
the canonical form/block-configuration associated with 
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FIG. 1: Evaluation of a 2-rdm element 74167. We can obtain 

this element e.g. at the block configuration where indices 4, 1 
are on the left block and indices 6, 7, are on the right block 
(corresponding to calling Compute(2, 0, 2) in Alg. 1). 



any site. By taking advantage of this flexibility, we 
can reduce the memory requirements once again back 
to 0{M'^k'^), i.e. the same as in the standard quantum 
chemical DMRG algorithm. Given a two-particle density 
matrix element {a''- (1^(1^(11) , where, say i < j < k < I, 
we choose a block configuration such that i,j lie in 
the left block and sites k,l lie in the right block, i.e. 

The correspond- 



ing matrix element may then be evaluated using ajaj^ on 



the left block, and 0^0; on the right block, and thus no 
operator matrices with more than two orbital indices ap- 
pear on either block (see Figure 1). By the appropriate 
choice of partitioning between the left and right blocks, 
wc can arrange things such that we never manipulate op- 
erators with more than two orbital labels on either the 
left or right blocks for any ijkl. During a DMRG sweep 
wc iterate through all block configurations where the di- 
viding site •p ranges from site 2 to site k — 1. At each 
block configuration, we then evaluate all the two-particle 
density matrix elements which do not require more than 
two- index operators on either the left or right blocks, and 
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assemble the contributions of all the block configurations 
at the end of the DMRG sweep. 

Along these lines, we can formulate an efficient algo- 
rithm to evaluate the two-particle density matrix with 
a total per-sweep computational cost of 0{M'^k^) and a 
memory cost of OiM^k"^). The pseudocode is given in 
Algs. (1), (2). Alg. (1) describes how to partition the 
evaluation of different density matrix elements amongst 
the block configurations as we traverse a DMRG sweep. 
The actual calculation of the density matrix elements is 
carried out by the function COMPUTE in Alg. (2), which 
computes all density matrix elements that may be assem- 
bled from nl index operators on the left block, np index 
operators on site p, and nr index operators on the right 
block. 

Algorithm 1 Two-particle density matrix evaluation 
showing how the two-particle density matrix is 
assembled across a DMRG sweep, 
special treatment for first configuration | •! | #2 



•3 . 



left= site 1, sitep= site 2, right= sites 3. 
Compute(4, 0, 0, left, siiep, right) 
Compute(3, 1, 0, left, sitep, right) 
Computers, 0, 1, left, sitep, right) 
Compute(2, 1, 1, left, sitep, right) 
sweep through block configurations 



•p+i . 



configuration 

right= site k 



for sitep— 2 to k-1 do 

left= sites 1 . . .p — 1, right= sites p + 1 
Compute(1, 2, 1, left, sitep, right) 
Compute(2, 1, 1, left, sitep, right) 
Compute(2, 2, 0, left, sitep, right) 
Compute(1, 3, 0, left, sitep, right) 
Compute(0, 3, 1, left, sitep, right) 
Compute(0, 4, 0, left, sitep, right) 

end for 

special treatment for final 

•l • • ■ •fc-2 'fc-l I 'fc I 

left= sites 1 ... fc — 2, sitep= site k - 
Compute(0, 0, 4, left, sitep, right) 
Compute(0, 1, 3, left, sitep, right) 
Compute(1, 0, 3, left, sitep, right) 
Compute(0, 2, 2, left, sitep, right) 
Compute(2, 0, 2, left, sitep, right) 
Compute(1, 1, 2, left, sitep, right) 
Compute(1, 2, 1, left, sitep, right) 



An attractive feature of the quantum chemical DMRG 
algorithm is the high level of parallelisability, which we 
have described in detail in Ref. [17]. In our imple- 
mentation, the loops over operators in Alg. (2) are 
trivially parallelised because of how our operators are 
divided across processors in our original formulation 
[17]. For example, the dominant computational cost of 
the two-particle density matrix evaluation comes from 
CoMPysTE{2,l,l, left, sitep, right) in Alg. (1), which 
costs 0{M^k^) per DMRG sweep. However, in our par- 
allel DMRG implementation, the two index operators 



Algorithm 2 COMP\JTE{nl,np,nr,left, sitep, right). 
Note nl, np, nr <2 and nl + np + nr = 4, i.e. the 
number of indices in the two-particle density matrix 7. 
for all opl= operators with nl indices on block left do 
(If parallel, loop only over opl stored on current proa) 
for all opp= operators with np indices on block sitep 
do 

for all opr= operators with nr indices on block right 
do 

7(np, nl, nr) — parity(op/, opp, opr) x {^\opl (g) opp® 
oprj^') 
end for 
end for 
end for 

(If parallel, accumulate contributions from all procs to root 
processor) 



opl on the left block, namely ajaj and ajflj, are divided 
across the processors, while the corresponding one in- 
dex operators opp, opr are replicated on all processors, 
and thus we can easily parallelise over the first opl loop 
in Alg. (2). This leads to a final computational cost 
per sweep of 0{M^k'^/np) with a communication cost of 
0(fc^ In rip), where rip is the number of processors. 



C. Orbital step and integral transformation 

As described earlier, the DMRG wavefunction is pri- 
marily efficient at capturing static correlation and conse- 
quently we employ an active space DMRG description of 
the electronic structure, the purpose of the orbital opti- 
misation then being to obtain the best form of the active 
space. Recall that the active space is defined by parti- 
tioning the orbitals into three sets, closed-shell orbitals 
which remain doubly occupied in all DMRG configura- 
tions, active orbitals which form the product active space 
{ni} (g) . . . (g) {rifc} in the DMRG wavefunction expansion 
(10), and external orbitals, which remain unoccupied in 
all DMRG configurations. With this partitioning, the 
active space DMRG wavefunction is determined with re- 
spect to the active space Hamiltonian 

^act ^ ^closed ^ ^ ^act^t^ . ^ ^ v.jkiaWjakai (25) 

ij ijkl 

where indices i,j are limited to the active orbitals and 

the modified one-particle integrals tfj^ and closed-shell 
energy are given respectively by 

^closed = J2tcc + $^(^^cc'c'c - t^cc'cc) (26) 

c cc' 

c 

where c, c' denote the closed-shell indices. 

Orbital optimisation chooses the best form of the ac- 
tive orbitals by minimising the energy of the DMRG 
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wavcfunction with respect to the active and closed-shell 
orbitals. This is the basic idea behind the Complete- 
Active-Space Self-Consistent Field (CASSCF) descrip- 
tion of electronic structure. In CASSCF, the active space 
wavefunction is the exact eigenfunction of the active 
space Hamiltonian (25) and is thus invariant with respect 
to active- active orbital rotations. In the corresponding 
orbital optimised DMRG-CASSCF, the accuracy of our 
active space DMRG wavefunction depends on the size of 
M, but in this study we will use sufSciently large M so 
that our wavefunction is nearly an exact eigenfunction of 
the active space Hamiltonian, and we will similarly omit 
active-active rotations. 

The algorithm we use for orbital optimisation is an 
Augmented Hessian Newton Raphson scheme similar to 
that used in modern CASSCF implementations [43, 44, 
50]. The orbital rotations are parameterised by the anti- 
hcrmitian amplitudes A in Eq. (5), and the derivative 
with respect to these amplitudes is evaluated from the 
one- and two-particle density matrices from the DMRG 
calculation. However, as the DMRG enables the use of 
larger active spaces than in traditional CASSCF studies 
and consequently we can expect to have a larger number 
of correlating external and closed-shell orbitals, we have 
focused on an efficient parallel implementation of the or- 
bital optimisation. Here the primary task is to paral- 
lelise the four-index transformation which is performed 
after each orbital rotation to generate the two-electron 
integrals in the basis of the rotated orbitals. We now 
describe how this is done. 

Say we have a coefficient matrix U giving the expan- 
sion coefficients for our rotated orbitals in terms of the 
starting atomic orbitals. Then, the transformed integrals 
Vpqrs axe obtained from the atomic orbital integrals 
through (assuming real coeflacients, for simplicity) 

VpQrs = Yl Up^,Ug,Ur^U,xV^Zx (28) 

As is well known, the four-index transformation should 

be carried out in four quarter-transformation steps corre- 
sponding to the four contractions with the coefficient ma- 
trices above. In our parallel transformation scheme, we 
consider the four steps in two stages; in the first stage we 
perform two quarter-transformations to construct half- 
transformed Coulomb and exchange intermediates J, K 

JaM'^) = Y.^a,U,xV^Zx (29) 
Kab{v, ^) = Y1 Ua^.UbXV^°^, (30) 

while in the second stage, we perform the remaining quar- 
ter transformations on the J, K intermediates to obtain 



the final integrals 

[Jabjpq = Vapqb = ^ Jabil^, K)Up„UqK, (31) 
[Kablpq = Vapbq = X! ^^"■bi'^^ K)Up„UgK (32) 

Note that for the purposes of optimising the active or- 
bitals, we only need the integrals that appear in the aug- 
mented Hessian. Thus, the ab indices in (29), (30) only 
need to run over the active orbitals while the pq indices 
need to run over all the closed-shell, active, and external 
orbitals. 

In the first stage, we parallelise the construction of 
the J, K intermediates by dividing up the intermediates 
according to their untransformed AO indices. For exam- 
ple, the construction of Jab{v, k) is divided amongst the 
processors according to the pair of indices {v, k) ; each 
processor is then responsible for constructing the J in- 
termediates for all {D, k) g proc. This allows us to also 
partition the AO integrals amongst the processors ac- 
cording to the same divided pair of indices {i^,K); e.g. 
to construct Jab{i^, «) for {v, k) G proc we only need AO 
integrals such as v^9.^ for {u, R) € proc to be stored on 
that processor. 

Once all J and K intermediates are constructed, we 
parallelise the second stage with respect to the trans- 
formed ab indices of the J, K intermediates. Thus ab is 
divided amongst the processors, and each processor con- 
structs the final integrals I'gpgj, I'apEg for all {ab} € proc. 
Since the first stage is parallelised over a pair of AO in- 
dices {v, k) (and the J and K intermediates are divided 
across the processors accordingly) while the second stage 
is parallelised over the two transformed indices {ab), we 
need to redistribute the intermediates J and K amongst 
the processors between the first and second stages. This 
is the main communication step. 

In addition to above parallelisation, further efficien- 
cies can be gained by using the permutational and spa- 
tial symmetries of the integrals. Our complete paral- 
lelised algorithm, which uses these symmetries, is pre- 
sented in pseudocode in Alg. (3). The cost of the 
four-index integral transformation as implemented is 
OiiK^k + K^k'^)/np) for CPU, OHK^ + K'^k'^)/np) for 
disk space, 0{K'^k'^ /up) for memory, and 0{K^k^) for 
overall communication, where K is the total number of 
orbitals, k is the number of active orbitals, and Up is the 
number of processors. 

To complete our efficient implementation of orbital 
optimisation, we have also parallelised the remaining 
steps in the Augmented Hessian Newton-Raphson solver. 
These additional steps take up only a small part of the 
computational time and have an overall cost 0{K^k^ /Up) 
for CPU time, 0{K^k^ /Up) for memory, 0{Kk) for com- 
munication. 
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Algorithm 3 Parallel four-index integral transformation 
algorithm. 

Stage 1: Assemble J and K intermediates 

Divide AO integrals by a factor (2-5^a)(2-(5,/«)(2- 

for i',k{i' > K) £ proc do 

for a, fx, X s.t. // > A, //A > t'/t do 

{i?, R) += Vf:g^ UaX ; m {i^, «) += <?^A Ua^. 

at; (17, n) += (7„A ; {y, k) += <° Ua^ 

end for 
for a, A do 

end for 

for a, /i, A s.t. fi > X, uii > fiX do 
end for 

for o, b, X s.t. o > 6 do 

Ll{v,R)Ut,x + L\{u,K)Uax 
end for 
for a, h, X do 

i^„(,(i>,A) +=iV^(i>,K)(76^ 
end for 
end for 

for a, h s.t. a > 6 do 

write Jab, Kab, and /sTba on disk 
end for 

Stage 2: Redistribute J and /sT, transform to final 

integrals 

for o, 6 (o > b) do 

road Jab, Kab, Kba from disk and send to proc(o, 6) 
end for 

for o, 6 (o > 6) e proc, z/, /t (f > /t) do 

Jab{l^,v)+ = Jab{y,R) 

end for 

for a,b{a >b) £ proc, i^, k do 
end for 

for a,b{a >b) £ proc, p, g, i^, k do 

f'spgE += Jab{'^,i^)Up^Ug^ (eqn. (31)) 

«dp6q += Kab{v,K)UpvUq^ (eqn. (32)) 
end for 



D. Complete Orbital Optimised DMRG-CASSCF 
Algorithm 

With the description of the density matrix evaluation 
in Sec. II B and the orbital optimisation and integral 
transformation in Sec. II C, we now have the basic in- 
gredients to perform the DMRG-CASSCF algorithm, ac- 
cording to the general outline in Sec. II A. 

There is one final ingredient however, the secret in- 
gredient. As the DMRG works best in a localised ba- 



sis (particularly in larger systems) it is beneficial to lo- 
calise the active space after each orbital optimisation. We 
have done this using the Pipck-Mezey procedure [51]; the 
active-space integrals are first transformed into this local 
basis before being input into the DMRG calculation. In 
total therefore, the complete DMRG-CASSCF algorithm 
is as follows: 

1. Localise the active space orbitals. 

2. Transform the AO integrals to the active space ba- 
sis and build the active space Hamiltonian. 

3. Perform the DMRG calculation using the active 
space Hamiltonian. 

4. From the converged DMRG wavefunctions at each 
block configuration, assemble the one- and two- 
particle density matrices. 

5. Using the density matrices, obtain the orbital gra- 
dient and orbital step from the Augmented Hessian 
Newton- Raphson solver. 

6. From the orbital step, determine the new active 
space orbitals. 

7. Goto 1. until convergence in the energy. 

Steps 1.-6. constitute a single DMRG-CASSCF macro- 
iteration. 



III. APPLICATIONS 
A. Long Polyenes 

1. Background 

Polyenes are the simplest conjugated systems, consist- 
ing of alternating singly and doubly bonded carbons ar- 
ranged in a chain. They are valuable models not only 
to understand conjugated polymers of materials inter- 
est (e.g. poly-acctylcne is simply an infinite polyene) 
but also biological molecules such as the carotenoid and 
retinal families of pigments involved in photosynthesis 
and vision. In these systems, the functionality of the 
molecules relics on the low-lying tt-tt* excited states of 
the conjugated backbone, which serve as the conduits 
for energy transfer. The excited states are labelled by 
their symmetry under the C2h point group, giving rise 
to Ag, Bg, Au, Bu symmetry labels. Furthermore, they 
are usually given an additional +/— label to indicate 
their approximate particle-hole symmetry. In Hamilto- 
nians (such as the Hiickel Hamiltonian) which support 
symmetric sets of energy states around the Fermi level, 
there is an additional symmetry associated with rotating 
the molecular orbital diagram so that the bonding and 
anti-bonding levels swap places [52]. Although particle- 
hole symmetry is not a true symmetry of the ab-initio 
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electronic Hamiltonian, it is still customary to use such 
labels for the polyenes, in particular, because the +/— 
states have very different qualitative electronic structure; 
valence bond studies of the Hubbard model [53] show 
that the + states consist mainly of ionic valence bond 
structures, while the — states consist mainly of covalent 
valence bond structures [42, 54, 55]. 

In this study we have looked only at singlet states and 
henceforth we shall be considering singlet states only. 
The ground state of the polyenes is known to always be of 
A~ symmetry. The lowest dipole-allowed singlet transi- 
tion, which has a predominantly HOMO^LUMO excita- 
tion character, has _B+ symmetry. However, contrary to 
what one might expect, this 1A~ 1_B+ transition is not 
the lowest singlet transition [56, 57]. Rather, as shown 
by Kohler et al. in octa-tetraene [56], there is a lower 
dipole forbidden excitation, later identified as the 2Ag 
state, which can be rationalised in valence bond language 
as arising from a pair of singlet-triplet excitations in the 
two separate double bonds that recouple to form a singlet 
state [58-66] . Following the observation of the 2A~ state 
in octa-tetraene, there has been much debate over the 
correct ordering of the 2A~ and 1-6+ excited states in 
the shorter polyenes, compounded both by experimental 
difficulties in observing the dipole-forbidden 2A~ state as 
well as theoretical challenges in achieving a balanced de- 
scription of the two states, which are dominated by very 
different kinds of correlation, namely static correlation in 
the 2A~ state and dynamic correlation in the 15+ state. 
In longer polyenes and the biologically active carotenoid 
and retinal pigments, questions about the low-lying spec- 
trum are not restricted simply to the 2A~ and 15+ state 
ordering. Recent studies using Resonance Raman exci- 
tation profiles (RREP) and electronic absorption spec- 
troscopy on substituted polyenes in the carotenoid fam- 
ily, have indicated the presence of additional dark states 
below the li?+ state [67 71]. In particular, for the all- 
iraKS-carotenoids with (the number of double bonds) 
n = 9 — 11, Sashima et al. observed a 1B~ state between 
the 2Ag and 1B+ [67, 72]. More recently, Furuichi et al. 
observed a 3^4^ level between the IB^ and 1_B+ states 
in carotenoids with n = 11 — 13, and assigned the tenta- 
tive state ordering of lA' < 2A- < IB' < SA' < 1B+ 
[71] . The assignment was made by extrapolating from the 
earlier PPP-MRDCI calculations by Tavan and Schulten 
on short polyenes (n = 2 — 8), which had predicted the 
existence of these additional states [55]. 

To better understand the electronic structure of these 
low-lying states, we would ideally like to be able to 
carry out an ab-initio multireference calculation, using 
the complete 7r-valence space. However, the large num- 
ber of active n orbitals in the longer polyenes means 
that it is not possible to perform such calculations with 
traditional CAS algorithms for these systems. Hirao 
and coworkers [41, 42] carried out incomplete valence 
CASSCF and CASCI-MRMP using a (10,10) active space 
on the polyene series up to C28H30 and observed rea- 
sonable agreement with experiment. However, with our 



new orbital optimised DMRG-CASSCF procedure, we 
can now re-examine the low-lying excitations in these 
systems correlating the complete tt- valence space even 
for the longer polyenes and carotenoids. 



2. Computational details 

The polyene molecular geometries for 

CsHio, C12H14, CigHis, C20H22, C24H26 were opti- 
mised at the density functional level using the B3LYP 
functional [73, 74] as implemented in Gaussian03 [75]. 
The polyene molecules were constrained to have C2h 
symmetry, with the C2 axis as the z-axis. The cc-pVDZ 
basis [76] was used for all calculations. 

In our DMRG-CASSCF calculations we used a com- 
plete TT- valence space i.e. in C24H26, this was a (24, 
24) active space. To generate this active space, wc first 
performed a restricted Hartree-Fock calculation in PSI3 
[77, 78] to obtain canonical Hartree-Fock molecular or- 
bitals. From these molecular orbitals, we could not triv- 
ially identify appropriate tt anti-bonding active orbitals 
because of significant 2p-3p mixing. We constructed the 
anti-bonding component of the active space as a set of 
projected atomic orbitals, by first projecting out the tt 
bonding space from a set of 2pz atomic orbitals. These 
projected atomic orbitals were then symmetrically or- 
thogonalised, then relocalised together with the bond- 
ing molecular orbitals (using the Pipek-Mczcy procedure 
[51]) to yield the complete active space in our calcula- 
tions. The final set of active orbitals generated in this 
way resemble an orthogonal set of 2pz orbitals. 

Note that our initial active space does not correspond 
precisely to an active space obtained by selecting Hartree- 
Fock canonical orbitals. Thus DMRG energies obtained 
before orbital optimisation do not correspond to typi- 
cal CASCI energies, but instead to CASCI energies ob- 
tained in our projected-atomic orbital (PAO) virtual 
space. This distinction is noted in our tables with the 
abbreviation DMRG-PAO-CASCI. After orbital optimi- 
sation, however, our DMRG-CASSCF energies do corre- 
spond to true CASSCF energies, up to the accuracy of 
the DMRG calculation. 

We carried out state-averaged DMRG-CASSCF calcu- 
lations in the above active space with the one-site DMRG 
algorithm with M = 250 and averaging over the 4 low- 
est eigenstates. The DMRG sweeps were converged to 
10~^"Eh in the DMRG energy, which took roughly 30 
DMRG sweeps. The number of renormalised states was 
increased smoothly from a starting value of M = 50 to 
the final value of M = 250. To aid the convergence of 
the DMRG sweeps in the one-site algorithm, we applied 
a system-environment perturbation as described in Ref. 
[79], with a starting magnitude of 10^'^ that smoothly 
decreased to after 20 sweeps. We estimate the remain- 
ing error in the DMRG energies at the M = 250 level 
from the exact Full-Configuration Interaction energies 
in the same active space to be less than Q.lmEh. Our 
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FIG. 2: DMRG-CASSCF excitation energies in eV for the 
2A~ , 1B~ and 3A~ states in the conjugated polyenes CsHio 

to C24H26 • 



FIG. 3: Converged DMRG sweep energies in Hartrees vs num- 
ber of orbital optimisation macro iterations in C20H22. 



DMRG calculations were combined with orbital rotation 
in a macro-iteration consisting of a converged DMRG cal- 
culation, an Augmented- Hessian step based orbital rota- 
tion, integral transformation, and orbital localisation, as 
described in Sec. II D. Typically 10-15 macro-iterations 
of the complete DMRG/orbital optimisation cycle were 
necessary to converge the energies to a tolerance of bet- 
ter than 10~^ Eh- The convergence of the state energies 
with the number of macro-iterations is shown in Fig. 3. 

The spatial and spin symmetries of excited states were 
assigned as follows. Firstly, all excited states were re- 
stricted to be of singlet spin symmetry through the ap- 
plication of a shift A(S'2 - (5) ((5) + 1) with A = 0.5 
[23]. To obtain the spatial symmetry, the ground state 
was assumed to be 1A~ as established by prior experi- 
mental and theoretical work. To determine whether the 
excited states were of Ag or symmetry the transi- 
tion dipole matrices were calculated between the states. 
Additionally, to determine the approximate particle-hole 
-|- or — symmetry we examined the magnitude of the 
transition dipoles; large transition dipoles for an allowed 
transition indicated that the transition involved a change 
of particle-hole symmetry between the states. 



3. Discussion 

In Table I we present the energies, symmetries, and 
oscillator strengths for the ground state and first 3 ex- 
citations in the polyenes from CsHio to C24H26- For 
comparison, we also give the excitation energies obtained 
from the CASCI-MRMP calculations of Kurashige et al. 
[42], as well as the experimental energies where avail- 




Size of Active Space 



FIG. 4: Change in CASSCF energies of the low-lying states 
of C12H14 as a function of increasing the active space from 
(4,4) to (12,12) (i.e. complete valence active space). 



able. (Note that in C20H22, the experimental excitation 
energies were obtained from the carotenoid spheroidene, 
which has a C20 conjugated backbone). 

We see that while our complete 7r-valence active space 
DMRG-CASSCF calculations generally overestimate the 
excitation energies, they reproduce the correct experi- 
mental ordering of the lowest excited states with the ex- 
ception of the missing IB+ state (the HOMO-LUMO ex- 
citation), which should lie helow the iA~ in the shorter 
polyenes such as CsHiq. If we perform a state-averaged 
DMRG-CASSCF with 5 states in CgHio, we find that 



11 



TABLE I: Energies, symmetries, and oscillator strengths for the lowest lying singlet excited states in conjugated polyenes. The 
DMRG-PAO-CASCI and DMRG-CASSCF entries for the lAj ground-states give the total energy in Eu; the other entries 
give the excitation energies from the ground state in eV. The estimated error of the DMRG-CASSCF energies from the exact 
CASSCF energies in the same active space is less than O.lmEh- The notation {n,m) denotes the active space used in the 
DMRG-PAO-CASCI and DMRG-CASSCF calculations. Oscillator strengths arc in a.u. for the ground-state, excited state 
transition. The CASCI-MRMP excitatiou energies are from Kurashigc et al. [42]; note that these used at most a (10,10) active 
space. The experimental numbers in brackets arc from measurements on the substituted polyene, spheroidene [71]. 
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the 1B+ state lies immediately above the 3A~. This 
may seem strange given that CASSCF is generally be- 
lieved to yield qualitatively correct electronic structure, 
but it reflects the wisdom from earlier studies on butadi- 
ene that fT-TT correlation is very strong in the 1J3+ state 
and must be included to obtain the correct balance be- 
tween Rydbcrg and valence character [65, 66, 82, 83]. 
Comparing with the calculations of Kurashige et al. [42] , 
which despite having an incomplete valence active space 
include dynamic ct-tt correlation through MRMP pertur- 
bation theory [84] , further indicates that ct-tt correlation 
would also lower the excitation energies of our other ex- 
cited states. 

To better understand the effect of using a complete tt 
valence space on the excitation energies, we have per- 
formed some small benchmark CASSCF calculations on 
C12H14 with 4 — 12 active orbitals. These results are 
prcscinted in Fig. 4. As can he seen, there; is a very 
strong dependence of the excitation energies on the size 
of the active space, and even the order of the excitations 
changes. Thus, while an incomplete valence active space 
can yield an excited state ordering in better agreement 
with experiment, one is tempted to argue that it does 
not do so for the right reason. 

In Fig. 5, we plot our DMRG-CASSCF excitation en- 



ergies as a function of the inverse chain length of the 
polyenes. Also shown (as an inset) is the same plot for 
the excitation energies obtained by Kurashige et al. [42] . 
It is easy to show that in a finite Hiickel model with n 
sites, the excitation energies have a sin(A;7r/2(2n -|- 1)) 
chain length dependence, where A; is a quasi-momentum 
number that labels the excitation. For long chains, this 
implies an asymptotic linear dependence on the inverse 
chain length 1 / (2n -|- 1) . Tavan and Schulten conjectured 
that this asymptotic behaviour held also in interacting 
systems, and presented evidence from MRD-CI calcula- 
tions on short-chain Hubbard (n up to 7) and Pariser- 
Parr-Pople models (n up to 8) to support the conjecture 
[85] . The experimental Resonance Raman excitation pro- 
flies from Sashima et al. [67] and Furuichi et al. [71] were 
also approximately fitted to the same inverse chain length 
behaviour, although only over a small range of n = 9—13. 
We see from our results that while the 2A~ and 1B~ ex- 
citation energies fit the asymptotic l/(2n-|- 1) behaviour 
well, the 3A~ state shows curvature more indicative of 
the sinusoidal dependence expected when k ^ 2n + 1. 
This is consistent with interpreting the 3A~ as an exci- 
tation labelled by a larger quasi-momentum than 2A~. 
Interestingly, the excitation energies of Kurashige et al. 
show quite different chain-length dependence, with all 
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FIG. 5: DMRG-CASSCF excitation energies for the low-lying 
singlet excited states of polyenes ranging from C12H14 to 

C24H26- The excitation energies are plotted against l/(2n-|~l) 
where n is the number of double bonds. The ratio of the 
slopes for the different states is found to be 2 : 3.0 : 3.8 as 
compared to 2 : 3.1 : 3.8 experimentally. Inset: same plot for 
the CASCI-MRMP excitation energies from Kurashige et al. 
[42]. As can be seen, these show a different and less linear- 
dependence on l/(2n -h 1). 



three states showing much stronger curvature when their 
excitation energies are plotted against l/{2n+ 1) in Fig. 
5 (inlay). Fitting our excitation energies for C16H20, 
C20H24, C24H26 (n = 8 — 12) to the asymptotic depen- 
dence l/(2n + 1), we obtain slopes of 27.67eV, 41.34eV, 
52.63eV for the 2A~, \B~ , iAg excitations, in reason- 
able agreement with the experimental slopes of 31.39eV, 
49.07eV and 59.63eV. 

From the one particle transition density matrices we 
can analyse the single-particle character of our exci- 
tations. Given the density matrix element Wij = 
(g.s.|a|aj|cxcitcd) where i, j arc natural orbitals in the 
ground state, we define the weight of the i ^ j excitation 
as w^j . The total single excitation weight is then ^ wfj . 
In Table II we give the largest excitation weights and the 
total single excitation weights for the low-lying polyene 
excited states as a function of the number of conjugated 
bonds. We see the 2A~, 1B~ and 3A~ states are domi- 
nated by many-particle excitations from the ground state 
(i.e. they have small single-particle excitation weights) 
and indeed the single-particle character of the excitations 
decreases even more as the chain-length increases. Re- 
markably, in C24H26 only < 16% of the excitation char- 
acter of these states can be considered to be of a single- 
particle nature! These results are consistent with the 
analysis by Wormer and Dreuw using coupled cluster and 
propagator techniques [86] . 
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TABLE II: Single particle nature of the polyene excitations 

(in %). For a given excited state (e.g. 2A~), the excitation 
weight of the transition i ^ j is given by [{lAg\alaj\2Ag)]^. 
The total excitation weight is the sum of weights for all transi- 
tions; 100% indicates that the given excited state corresponds 
entirely to single excitations from the ground state. The tran- 
sition labels n —urn' are interpreted as follows: 1, 2, 3 ... de- 
note HOMO, HOMO-1, HOMO-2 . . . natural orbitals, while 
1', 2', 3' denote LUMO, LUMO+l,LUMO+2 natural orbitals. 
As the polyenes increase in length, the total weight of the 
single excitations in the low-lying states becomes very small, 
< 16%. 



State 


Excitation 


No. of conjugated double bonds 






weight 


4 


6 


8 


10 


12 


2A; 


2^1' 


10.9 


8.6 


6.6 


5.3 


4.3 




1-^2' 


6.7 


5.9 


4.8 


4.0 


3.3 




Total 


20.0 


18.0 


15.4 


13.5 


12.1 




3^1' 


14.5 


10.2 


7.9 


6.3 


5.2 




1^3' 


7.0 


5.6 


4.6 


3.9 


3.3 




Total 


25.3 


21.8 


18.6 


16.3 


14.7 


3A- 


4^1' 


21.3 


12.8 


9.3 


7.1 


5.6 




1^4' 


8.2 


6.0 


4.7 


3.8 


3.1 




Total 


32.9 


25.0 


20.9 


18.0 


15.9 



B. /3-carotene 




FIG. 6: s-cis /3-carotene. 



Carotenoids, the family of substituted polyenes, are 
the primary light harvesting pigments in the LH2 com- 
plex. Light harvesting proceeds by the transfer of energy 
from an array of carotenoids to nearby bacteriochloro- 
phylls and thence to the photosynthctic centre. Many 
essential questions remain unanswered as to the precise 
mechanism of this energy transfer [86-92]. While the 
absorption of light places the carotenoid in the dipole 
allowed excited state, there can be a fast internal con- 
version to the aforementioned dark states of the polyene 
backbone, and thus multiple pathways for energy transfer 
to the bacteriochlorophyll. In carotenoids, the dipole al- 
lowed transition is usually labelled 5*2, while historically 
the dark state is labelled 5*1. However, with the dis- 
covery, as previously described, of additional dark states 
below 5*2 in these molecules [67-72], this nomenclature 
can be confusing. An alternative nomenclature is to sim- 
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7.0 




2.5-1 1 . 1 . 1 • 1 • — 

4 6 8 10 12 

No. of conjugated double bonds 



FIG. 7: Polyene and carotene excitation energies vs the num- 
ber of double bonds: the /3-carotene excitation energies when 
fitted to the polyene excitation energies give an effective con- 
jugation length of 9.5 — 9.7. 



ply re-use the polyene excited state labels, even though 
the carotenoids have a lower point group symmetry We 
will follow this practice here. 



1. Discussion 



TABLE III: DMRG-CASSCF energies, symmetries, and os- 
cillator strengths for the lowest lying singlet excited states in 
/3-carotene with the complete 7r-valence (22,22) active space. 
Total energies in Eh, excitation energies in eV, oscillator 
strengths in a.u.. The estimated error of the DMRG-CASSCF 
energies from the exact CASSCF energies in the same active 
space is less than O.lmEh- Oscillator strengths are for the 
ground-state, excited state transition. 



Symmetry 


DMRG-CASSCF 


Excitation 


Oscillator 


Expt 




total energy 


energy 


Strength 




lA- 


-1546.914545 








2A- 


-1546.804503 


2.99 


Forbidden 


1.81 " 


IB- 


-1546.781125 


3.63 


0.2025 


2.05 " 


3A- 


-1546.755822 


4.31 


Forbidden 


(2.22) ' 



"[68]. 

''Excitation measured for lycopene [71]. 



We have chosen to study s-cis /3-carotene (see Fig. 6) 
as a representative carotenoid. It is the dominant natu- 
ral conformer although the all-trans form is also studied. 
Crystalline /3-carotcne has Ci symmetry with a conju- 
gated backbone that lies almost entirely on the xy plane 
except for end groups which are twisted out of plane 




(a)LUMO-l-l natural orbital 




(b)LUMO natural orbital 




(c)HOMO natural orbital 



.J? 




(d)HOMO-l natural orbital 

FIG. 8: Natural orbitals corresponding to the HOMO-1 
through LUMO-l-1 states. 

These orbitals participate in the lowest lying singlet exci- 
tations in /3-carotene and contain little density on the non- 
planar end groups. 



[93, 94]. (In the biological setting, carotenoid pigments 
usually adopt a twisted configuration in the conjugated 
backbone[95, 96]). There are 11 conjugated double bonds 
in the backbone. Our study employed the same calcu- 
lation procedure as described in Sec. Ill A 2 with the 
exception that we used a 6-31G basis set in the DMRG- 
CASSCF calculation due to the large size of the molecule. 
State-averaged DJVIRG-CASSCF calculations were per- 
formed with 4 states and a (22,22) complete 7r-valence 
space, in the manner described in Sec. Ill A 2. 

In Table III we present the energies, symmetries, and 
oscillator strengths for the ground state and first 3 exci- 
tations in /3-carotene. We reproduce the state ordering 
lAg < 2Ag < li?,j7 < 3A~ as assigned by Furuichi et 
al. [71] (note that the 1B+ which does not appear in 
our calculation indeed lies above the 3^4" state in this 
molecule). However, just as in the polyenes, the exci- 
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tation energies from the DMRG-CASSCF procedure are 
generally overestimated in comparison with experiment, 
most likely due to the lack of a--K dynamic correlation. 

A question that has received some attention in the lit- 
erature is the effective conjugation length of carotenoids, 
since the presence of substituents and non-planar geome- 
tries are expected to modify this from the naive value 
deduced from the Lewis structure [97]. Formally, /3- 
carotene has 11 double bonds in the polyene backbone, 
but by comparing the excitation energies of the polyenes 
with our /3-carotene excitation energies, we can estimate 
a reduced conjugation length of 9.5-9.7 bonds, which is 
very close to the experimental estimate of 9.7 of Onaka 
et al. [70] . This reduced conjugation length results from 
the twist in the carotene end-groups. In Fig. 8 we plot 
the DMRG-CASSCF natural orbitals corresponding to 
the HOMO, HOMO-1, LUMO, and LUMO-fl. As can 
be seen, there is very little density in these orbitals on 
the carotene end-groups, and this is consistent with our 
reduced effective conjugation length. 

IV. CONCLUSION 

In this work, we described how to efficiently implement 

orbital optimisation using the Density Matrix Renormal- 
ization Group (DMRG) wavefunction. We have named 
the resulting method DMRG-CASSCF, and by virtue of 
the compact nature of the DMRG wavefunction, this now 
enables us to handle much larger active spaces than are 



possible with the traditional CASSCF algorithm. As a 
sample application, we have used our DMRG-CASSCF 
implementation to study the low-lying excitations of 
polyenes from CsHio to C24H26 as well as the light- 
harvesting pigment /3-carotene, with up to a (24,24) com- 
plete active space. Our calculations reproduce the state 
ordering of the dark states that have been recently ob- 
served by Resonance Raman studies. However, as ex- 
pected from earlier CASSCF studies, the energy of the 
optically allowed HOMO-LUMO 1B+ transition is still 
overestimated, as a result of the lack of dynamic ct-tt 
correlation in the DMRG-CASSCF method. We there- 
fore view the incorporation of dynamic correlation, either 
via perturbation theory or via canonical transformation 
[98, 99] into the DMRG-CASSCF method to present an 
important next direction for development. 
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